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Abstract. Numerically solving the Boltzmann kinetic equations with the small Knudsen number 
is challenging due to the stiff nonlinear collision term. A class of asymptotic preserving schemes 
was introduced in [6] to handle this kind of problems. The idea is to penalize the stiff collision 
term by a BGK type operator. This method, however, encounters its own difficulty when applied 
to the quantum Boltzmann equation. To define the quantum Maxwellian (Bose-Einstein or Fermi- 
Dirac distribution) at each time step and every mesh point, one has to invert a nonlinear equation 
that connects the macroscopic quantity fugacity with density and internal energy. Setting a good 
initial guess for the iterative method is troublesome in most cases because of the complexity of the 
quantum functions (Bose-Einstein or Fermi-Dirac function). In this paper, we propose to penalize 
the quantum collision term by a 'classical' BGK operator instead of the quantum one. This is based 
on the observation that the classical Maxwellian, with the temperature replaced by the internal 
energy, has the same first five moments as the quantum Maxwellian. The scheme so designed avoids 
the aforementioned difficulty, and one can show that the density distribution is still driven toward 
the quantum equilibrium. Numerical results are present to illustrate the efficiency of the new scheme 
in both the hydrodynamic and kinetic regimes. We also develop a spectral method for the quantum 
collision operator. 

1. Introduction 

The quantum Boltzmann equation (QBE), also known as the Uehling-Uhlenbeck equation, describes 
the behaviors of a dilute quantum gas. It was first formulated by Nordheim [T3] and Uehling and 
Uhlenbeck [16 from the classical Boltzmann equation by heuristic arguments. Here we mainly consider 
two kinds of quantum gases: the Bose gas and the Fermi gas. The Bose gas is composed of Bosons, 
which have an integer value of spin, and obey the Bose-Einstein statistics. The Fermi gas is composed 
of Fermions, which have half-integer spins and obey the Fermi-Dirac statistics. 

Let f(t,x,v) > be the phase space distribution function depending on time t, position x and 
particle velocity v, then the quantum Boltzmann equation reads: 

(1.1) y-+v-V x f=^Q q (f), leSJc! 4 ,^! 4 '. 

Here e is the Knudsen number which measures the degree of rarefaction of a gas. It is the ratio 
between the mean free path and the typical length scale. The quantum collision operator Q q is 

(1.2) Q q (f)(v)= f I B{v-v*,oj)[ffi{i±e f){i±e Q u)-fU{i±e f'){i±e fi)]du>dv* 

where 8q = h dv , H is the rescaled Planck constant. In this paper, the upper sign will always correspond 
to the Bose gas while the lower sign to the Fermi gas. For the Fermi gas, we also need / < by 
the Pauli exclusion principle. /, /*, /' and fl are the shorthand notations for f(t,x,v), f(t,x,v*), 
f(t,x,v') and f(t,x,v^) respectively, (v, «*) and (i>',i>*) are the velocities before and after collision. 
They are related by the following parametrization: 

, v + v* \v — v*\ 



(1.3) 



v + v* \v — v*\ 

7. ^ w ) 
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where u> is the unit vector along v' — v' t . The collision kernel B is a nonnegative function that only 
depends on \v — v*\ and cos 9 (9 is the angle between ui and v — u*). In the Variable Hard Sphere 
(VHS) model, it is given by 

(1.4) B(v-v*,u) = C 1 \v-v*\" 1 

where C 7 is a positive constant. 7 = corresponds to the Maxwellian molecules, 7 = 1 is the hard 
sphere model. 

When the Knudsen number e is small, the right hand side of equation (jl.ip becomes stiff and 
explicit schemes are subject to severe stability constraints. Implicit schemes allow larger time step, 
but new difficulty arises in seeking the numerical solution of a fully nonlinear problem at each time 
step. Ideally, one wants an implicit scheme allowing large time steps and can be inverted easily. In 
[5], for the classical Boltzmann equation, Filbet and Jin proposed to penalize the nonlinear collision 
operator Q c by a BGK operator: 

(1.5) Q c = [Q c - X(M C - /)] + \[M C - /] 

where A is a constant that depends on the spectral radius of the linearized collision operator of Q c 
around the local (classical) Maxwellian M. c . Now the term in the first bracket of the right hand side 
of (jl.5[) is less stiff than the second one and can be treated explicitly. The term in the second bracket 
will be discretized implicitly. Using the conservation property of the BGK operator, this implicit 
term can actually be solved explicitly. Thus they arrive at a scheme which is uniformly stable in e, 
with an implicit source term that can be inverted explicitly. Furthermore, under certain conditions, 
one could show that this type of schemes has the following property: the distance between / and 
the Maxwellian will be 0(e) after several time steps, no matter what the initial condition is. This 
guarantees the capturing of the fluid dynamic limit even if the time step is larger than the mean free 
time. 

Back to the quantum Boltzmann equation (II. lj) . a natural way to generalize the above idea is to 
penalize Q q with the quantum BGK operator Ai q — /. This means we have to invert a nonlinear 
algebraic system that contains the unknown quantum Maxwellian A4 q (Bose-Einstein or Fermi-Dirac 
distribution) for every time step. As mentioned in [7J, this is not a trivial task compared to the 
classical case. Specifically, one has to invert a nonlinear 2 by 2 system (can be reduced to one nonlinear 
equation) to obtain the macroscopic quantities, temperature and fugacity. Due to the complexity of 
the quantum distribution functions (Bose-Einstein or Fermi-Dirac function) , it is really a delicate issue 
to set a good initial guess for an iterative method such as the Newton method to converge. 

In this work we propose a new scheme for the quantum Boltzmann equation. Our idea is based 
on the observation that the classical Maxwellian, with the temperature replaced by the (quantum) 
internal energy, has the same first five moments as the quantum Maxwellian. This observation was 
used in [JJ to derive a 'classical' kinetic scheme for the quantum hydrodynamical equations. Therefore, 
we just penalize the quantum collision operator Q q by a 'classical' BGK operator, thus avoid the 
aforementioned difficulty. At the same time, we have to sacrifice a little bit on the asymptotic property. 
Later we will prove that for the quantum BGK equation, the so obtained / satisfies: 

(1.6) /" - M q = O(At) for some n > N, any initial data f°, 

i.e. / will converge to the quantum Maxwellian beyond the initial layer with an error of O(At). 

Another numerical issue is how to evaluate the quantum collision operator Q q . In fact (ll.2[) can be 
simplified as 

(1.7) Q q (f)(v)= f f B(v-v*,u)[f , fUl±e f±6 f*)-ff*(l±6 f'±d fi)]<h]dv* 

so Q q is indeed a cubic operator. Almost all the existing fast algorithms are designed for the classical 
Boltzmann operator based on its quadratic structure. Here we will give a spectral method for the 
approximation of Q q . As far as we know, this is the first time to compute the full quantum Boltzmann 
collision operator with the spectral accuracy. 

The rest of the paper is organized as follows. In the next section, we give a brief introduction to the 
quantum Boltzmann equation: the basic properties, the quantum Maxwellians and the hydrodynamic 
limits. In section 3, we present the details of computing the quantum collision operator by the spectral 
method as well as the numerical accuracy. Our new scheme to capture the hydrodynamic regime is 
given in section 4. In section 5, the proposed schemes are tested on the 1-D shock tube problem of 
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the quantum gas for different Knudsen number e ranging from fluid regime to kinetic regime. The 
behaviors of the Bose gas and the Fermi gas in both the classical regime and quantum regime are 
included. Finally some concluding remarks are given in section 6. 



2. The Quantum Boltzmann Equation and its Hydrodynamic Limits 
In this section we review some basic facts about the quantum Boltzmann equation (jl.ljl . 
• At the formal level, Q q conserves mass, momentum and energy. 



(2.1) 



/ Q q (f)dv= [ Q q (f)vdv= f Q q (f)\v\ 2 dv = 0. 

JR d v jRd v J R dv 

If / is a solution of QBE (jl.ip . the following local conservation laws hold: 
d 



(2.2) 



oil fdv + v > 



d 



d_ 

I at 



vfdv + V x 
l\v\ 2 fdv + \7 :i 



vfdv = 0, 
v <8> vfdv = 0, 

1 



v—\v\ 2 fdv = 0. 



(2.3) 



(2.4) 



Define the macroscopic quantities: density p, macroscopic velocity u, specific internal energy 
e as 

P= [ fdv, pu= I vfdv, pe= I ]-\v-u\ 2 fdv 

Jg.d v J R d v J M d v I 

and stress tensor P and heat flux q 

(v — u) ® (v — u)fdv, q 

the above system can then be recast as 
( dp 



-{v - u)\v - u\ 2 fdv, 



(2.5) 



at 

d(pu) 
dt 



+ V x ■ (pu) = 0, 
+ \7 X ■ (pu ® u + '. 



• Q q satisfies Boltzmann's H-Theorem, 



(2.6) 



(2.7) 



(2.8) 



In 



/ 



l±0o/ 



pe + -pu z ) u + Pit + q ) = 0. 



Q q (f)dv < 0, 



moreover, 



In 



/ 



l±0of 



Q q {f)dv = ^ Q q (f) =0^f = M q , 



where Ai q is the quantum Maxwellian given by 



M 



1 



1 



z _1 e 2T zp 1 



where z is the fugacity, T is the temperature (see [7] for more details about the derivation of 
A4 q ). This is the well-known Bose-Einstein ('-') and Fermi-Dirac ('+') distributions. 
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2.1. Hydrodynamic Limits. Substituting M. q into 
yielding the quantum Euler equations: 



^ + V x ■ (pu) = 0, 



(2.9) 



d(pu) 
dt 



V x ■ [ pu ® u + -j-P e I 



0. 



3^ 



2 PU 



-pe 



(|2.4I) . the system (|2.5|) can be closed, 



2 PU 



u = 0. 



With the macroscopic variables p, u and e, they are exactly the same as the classical Euler equations. 
However, the intrinsic constitutive relation is quite different, p and e are connected with T and z 
(used in the definition of M q (|2.8p ) by a nonlinear 2 by 2 system: 



(2.10) 



2 Q^(z) ' 



where Q u (z) denotes the Bose- Einstein function G v (z) and the Fermi-Dirac function -Fly(z) respec- 
tively, 



(2.11) 
(2.12) 



G„(z) 



1 



i»./ z-^e* 



1 



r(^) y z-^ + i 



dx, < z < 1, f > 0; 2 = 1, u > 1, 



dx, < z < oo, ^ > 0, 



and r(^) = Jq 00 x v ~ 1 e~ x dx is the Gamma function. 

The physical range of interest for a Bose gas is < z < 1, where z = 1 corresponds to the degenerate 
case (the onset of Bose-Einstein condensation). For the Fermi gas we don't have such a restriction 
and the degenerate case is reached when z is very large. For small z (0 < z < 1), the integrand in 
(|2.1ip and (|2.12[) can be expanded in powers of z, 



(2.13) 
(2.14) 



n=l 



z 2 z 3 



OO 

E 



(-l)" +i — 



2" 



z 



Thus, for z < 1, both functions behave like z itself and one recovers the classical limit. 
On the other hand, the first equation of (|2.10l) can be written as 



(2.15) 



.(*) 



is just the coefficient of the classical Maxwellian, which should be an O(l) quantity. 



where ^-3— 

(2ttT)^ L 

Now if 6*o —5- 0, then Qd 3L (z) — > 0, which means z <C 1 by the monotonicity of the function Q v . This 
is consistent with the fact that one gets the classical Boltzmann equation in QBE (|1.1[) by letting 
O -> 0. 

The quantum Euler equations (|2.9p can be derived via the Chapman-Enskog expansion [3j as the 
leading order approximation of the quantum Boltzmann equation (11.11) . By going to the next order, 
one can also obtain the quantum Navier-Stokes system which differs from their classical counterparts. 
In particular, the viscosity coefficient and the heat conductivity depend upon both p and e pQ. 
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3. Computing the Quantum Collision Operator Q q 

In this section, we discuss the approximation of the quantum collision operator Q q . The method 
we use is an extension of the spectral method introduced in [12l [5] for the classical collision operator. 
We first write (|1.2[) as 

(3.1) Q q = Qc±0o{Qi + Q2-Q3-Qi), 

where 

(3.2) 

is the classical collision operator. The cubic terms Qi - Q4 are 



Qc(f )(v) = I / B(v- v„u)\ff m - ff*]dojdv 



(3-3) 



&(/)(«) = 
Ca (/)(«) = 
Ca(/)(«) = 



§d v -l 



B(v - v*,u))f fj*dudv, 



B(v - v*,ui)f'flfdwdv, 



B(v - v*,uj)ff*f'dudv, 



B(v - v*,u))ff*fldujdv. 



In order to perform the Fourier transform, we periodize the function / on the domain T>l = 
[-L, L] dv (L is chosen such that L > 3+ 2 v ^ R, R is the truncation of the collision integral which 
satisfies R = 2S, where B(Q,S) is an approximation of the support of / [Hj). Using the Carleman 
representation [2] , one can rewrite the operators as (for simplicity we only consider the 2-D Maxwellian 
molecules), 



(3.4) 
and 



Qc(f)(v)= / / S(x-y)[f{v + x)f(v + y)-f(v + x + y)f(v)]dxdy 

J Br J Br 



(3.5) 



&(/)(«) = 
&(/)(«) = 
&(/)(«) = 
Q 4 (/)(«) = 



Bh 



Br 



<5(^ ' y)fi v + x)f(v + y)f(v + x + y)dxdy, 



K x ■ y)f( v + x )f( v + y)f(v)dxdy, 



&{ x ' V)f( v + x )f( v + x + y)f(v)dxdy, 



H x ■ y)f{v + y)f(v + x + y)f(v)dxdy. 



'Br JBr 

Now we approximate / by a truncated Fourier series, 
(3.6) /(«)« J2 heW, h 



1 



fe=- 



(2L) 



f{v)e-^""dv 



Plugging it into (|3.4|l (|3.5|l . one can get the fc-th mode of Q q . The classical part is the same as those 
in the previous method |12j . We will mainly focus on the cubic terms. 
Define the kernel modes 



(3.7) 



/3(/,m) 



Br JB 



5(x ■ y)e t i lx e l i m - y dxdy. 
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Following 

(3.8) 

with 
(3.9) 

where d>(s) 



], j3(l, m) can be decomposed as 



M-l 

P(l,m) = — E OL p {l)a' p (m) 



J9=0 



2L 

7TS 



a p (l) — 4>{l ■ (cos 6 p , sm.6 p )), ct' p (m) = <f>(m ■ (— sin(9 p , cos^p)), 
sin( fj-Rs), M is the number of equally spaced points in [0, and 6 P = f 



Then 



• The fe-th coefficient of Qi is 

2 1 

^ + n,m + n)fif m f n = 

l,m,n——M- 
l-\-m-\-n—k 



M-l f-1 

-Y Y 



M 



p=0 



Y. a p (l + n)a' p (m + n)fif„ 



i,m=-f 
l+7n—k—n 



M-l 



(3.10) 



E 9k-n{n)f„ 



M 



Terms inside the bracket is a convolution (defined as gk-n(n)), which can be computed by 
the Fast Fourier Transform (FFT). However, the outside structure is not a convolution, since 
9k-n{n) itself depends on n. So we compute this part directly. 
• The k-th coefficient of Q2 is 



M-l 



(3.11) 



l,m,Tl--y 

l+m+n—k 



p=0 



E a p( l ) a p( m )flfrr 



1+m—k—n 



fn- 



In this case, both inside and outside are convolutions. The FFT can be implemented easily. 
• The k-th coefficient of Q3 is 



M-l 



(3.12) e w+'^h/iu^E E a p( l 



1 iV 

L,m,n— — 
Z+m+n— k 



p=0 



a p( m )flfn 



l-^rn—k—n 

Factoring out a p (l 4- m), both inside and outside are convolutions again. 
The k-th coefficient of 0,4 is 



fn • 



E a p( m )fi.u 



l+m—k—n 



f-1 M-l f-1 

(3.13) E P{m,l + m)fifJ n = ^Y. E ^ + ' 

i,m,n=-f P=0 „=-f 

/+m+n— /e 

This term can be evaluated similarly as Q3. 

Remark 3.1. The computational cost of this quantum solver is 0(MN A \ogN), which mainly comes 
from computing Q\. This cost is higher than 0(MN A ) of the discrete velocity model. But taking 
into account the high accuracy and small value of log N (N is not very big in the real simulation), 
our method is still more attractive than the quadrature method. The fast algorithm for the quantum 
collision operator remains an open problem. 

3.1. Numerical Accuracy. To illustrate the accuracy of the above method, we test it on a steady 
state, namely, we compute Q q (M. q ) and check its max norm. In all the numerical simulations, the 
particles are assumed to be the 2-D Maxwellian molecules. 

Let p = 1, T = 1, from (|2.10p one can adjust 6q to get z that lies in different physical regimes. 
When 9 = 0.01 (h = 0.1), zbosc = 0.001590, ZFermi = 0.001593. In this situation, the quantum effect 
is very small. The Maxwellians for the Bose gas, classical gas and Fermi gas are almost the same 
(FigfTJ). When we increase #0, sa Y #0 = 9 (h — 3), zbosc = 0.761263, ZFcrmi = 3.188717, the difference 
between the quantum gases and the classical gas is evident (Figf2J). 
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Bose classical Fermi 




Figure 1. The Maxwellians at p = 1, T = 1, Oq = 0.01. Left: Bose gas; Center: 
classical gas; Right: Fermi gas. 



Bose classical Fermi 




Figure 2. The Maxwellians at p = 1, T = 1, 8 = 9. Left: Bose gas; Center: 
classical gas (same as in FigfT]); Right: Fermi gas. 

In Tabled! we list the values of || Q c {M-c) \\l°° and || Q q {M. q ) \\l<*> computed on different meshes 
N=16, 32, 64 (number of points in v direction), M=4 (number of points in angular direction 9 p ; it 
is not necessary to put too many points since M won't effect the spectral accuracy, see [H]). The 
computational domain is [—8,8] x [—8,8] (L = 8). 

These results confirm the spectral accuracy of the method, although the accuracy in the quantum 
regime is not as good as that in the classical regime. This is because the regularity of the quantum 
Maxwellians becomes worse when 6q is increasing, or strictly speaking, the mesh size Av is not small 
enough to capture the shape of the Maxwellians. To remedy this problem, one can add more grid 
points or more effectively, shorten the computational domain. For the Bose-Einstein distribution, we 
also include the results computed on [—6, 6] x [—6, 6] in TableJTJ One can clearly see the improvements. 
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16 x 16 


32 x 32 


64 x 64 


convergence rate 


classical gas 


2.1746e-04 


3.8063e-12 


1.9095e-16 


20.0253 


Bose gas 9 = 0.01 
0o = 9 

O = 9, L = 6 


2.1084e-04 
0.4891 
0.1815 


2.5512e-10 
0.0310 
0.0052 


1.9080e-16 
1.3496e-04 
4.0278e-06 


20.0036 
5.9117 
7.7298 


Fermi gas do = 0.01 
0o = 9 


2.2397e-04 
8.9338e-04 


1.6485e-10 
2.0192e-06 


1.9152e-16 
1.5962e-10 


20.0445 
11.2081 



Table 1. Comparison of the quantum collision solver on different Maxwellians (L = 8 
unless specified). 



3.2. Relaxation to Equilibrium. Let us consider the space homogeneous quantum Boltzmann equa- 
tion for the 2-D Maxwellian molecules. As already mentioned, this equation satisfies the entropy con- 
dition, and the equilibrium states are the entropy minimizers. Hence, we first consider the quantum 
Boltzmann equation for a Fermi gas with an initial datum < fo < and observe the relaxation 
to equilibrium of the distribution function. Then, we take a Bose gas for which the entropy is now 
sublinear and fails to prevent concentration, which is consistent with the fact that condensation may 
occur in the long-time limit. 

Fermi gas. The initial data is chosen as the sum of two Maxwellian functions 



(3.14) 



fo(v) = exp 



exp 



- Vi\ 



v e 



with v\ = (2, 1). The final time of the simulation is T en d = 0.5, which is very close to the stationary 
state. 

In the spatially homogeneous setting, Pauli's exclusion principle facilitates things because of the 
additional L°° bound < f(t) < j-. In this case, the convergence to equilibrium in a weak sense 
has been shown by Lu 10 j. Later Lu and Wennberg proved the strong L 1 stability [5]. However, 
no constructive result in this direction has ever been obtained, neither has any entropy-dissipation 
inequality been established. 

In Fig(3] we report the time evolution of the entropy and the fourth and sixth order moments of 
the distribution with respect to the velocity variable. We indeed observe the convergence to a steady 
state of the entropy and also of high order moments when t — > oo. 




1.25 



fourth order moment 



sixlh order mrmEnt 



0.1 0.2 0.3 OA 0.5 




C 0.1 2 0.3 0.4 OA 




Figure 3. Fermi gas. Time evolution of the entropy, fourth and sixth order moments. 

In Fig|4]we also report the time evolution of the level set of the distribution function f(t,v x ,v y ) 
obtained with N = 64 modes at different times. Initially the level set of the initial data corresponds 
to two spheres in the velocity space. Then, the two distributions start to mix together until the 
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stationary state is reached, represented by a single centered sphere. It is clear that the spherical 
shapes of the level sets are described with great accuracy by the spectral method. 



t=0.00 t=0.02 




Figure 4. Fermi gas. Time evolution of the distribution function f(t,v x ,v y ) with 
N = 64 modes at times t = 0, 0.02, 0.04 and 0.5. 

Bose gas. This is an even more challenging problem since there is no convergence result, due to the 
lack of a priori bound. Lu [TT] has attacked this problem with the well-developed tools of the modern 
spatially homogeneous theory and proved that the solution (with a very low temperature) converges 
to equilibrium in a weak sense. In [4], the authors studied an one dimensional model and proved 
existence theorems, and convergence to a Bose distribution having a singularity when time goes to 
infinity because Bose condensation cannot occur in finite time. 

Here we investigate the convergence to equilibrium for space homogeneous model in 2-D, for which 
condensation cannot occur. We consider the following initial datum 

(3.15) W ,) = ^^-t^)+» P (-^) ; .etf, 

with vi = (1, 1/2) and T = 1/4. 

We still observe the convergence to equilibrium and convergence of high order moments when t — > oo 
in Fig|5] 

In Figj6]we report the time evolution of the level set of the distribution function f(t, v x ,v y ) obtained 
with N = 64 modes at different times and observe the trend to equilibrium. 

4. A Scheme Efficient in the Fluid Regime 

So far we have only considered spatially homogeneous quantum Boltzmann equations, now what 
happens for spatially inhomogeneous data? Due to the natural bound < f(t) < 4-, the Boltzmann- 
Fermi model seems to be well understood mathematically |17j . The situation is completely different 
for the Boltzmann-Bose model, since singular measures may occur |17j . 
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fin-ropy 



lourih under momm 



sixth ord=r moment 





Figure 5. Bose gas. Time evolution of the entropy, fourth and sixth order moments. 



tO .00 MMK 




Figure 6. Bose gas. Time evolution of the distribution function f(t,v x ,v y ) with 
N = 64 modes at times t = 0, 0.02, 0.04 and 0.5. 

We first review the scheme in [5] for the classical Boltzmann equation 

(4.1) % +v-V x f=-Q c (f). 

at e 

The first-order scheme reads: 

(4.2) -—— +v . Vx f = . + . , 

where A is some appropriate approximation of |VQ C | (can be made time dependent). To solve f n+1 
explicitly, we need to compute Ai" +1 first. Since the right hand side of (|4.2p is conservative, it vanishes 



A NUMERICAL SCHEME FOR THE QUANTUM BOLTZMANN EQUATION EFFICIENT IN THE FLUID REGIME 1 



when we take the moments (multiply by <fi(v) — (l,v, \v 2 ) T and integrate with respect to v). Then 
(|4.2D becomes 

(4.3) — + / 4>(v)v ■ V x f n dv = 0, 

where U = (p, pu, pe + ^pu 2 ) T is the conserved quantities. Once we get U n+1 , M™ +1 is known. Now 
f n+1 in (l4~2l) is easy to obtain. 

When generalizing the above idea to the quantum Boltzmann equation (jl.ip , the natural idea is to 
replace Q c and M. c in (|4.2j) by Q q and M. q respectively. However, as mentioned in section 2, one has 
to invert the nonlinear system (|2.10j) to get z and T. Experiments show that the iterative methods 
do converge when the initial guess is close to the solution (analytically, this system has a solution [T]). 
But how to set a good initial guess for every spatial point and every time step is not an easy task, 
especially when p and e are not continuous. 

Here we propose to use a 'classical' BGK operator to penalize Q q . Specifically, we replace the 
temperature T with the internal energy e in the classical Maxwellian using relation e = %-T (true for 
classical monatomic gases) and get 

An important property of A4 C is that it has the same first five moments as Ai q . 
Now our new scheme for QBE (jl.ip can be written as 

(4 5) - /" , v Vxfn _ QM n ) - KM n c - r) | \{M n c +1 - /"+ 1 ) _ 

At x e e 

Since the right hand side is still conservative, one computes M™ +1 the same as for (|4.2j) . 

It is important to notice that z and T are not present at all in this new scheme, thus one does not 
need to invert the 2 by 2 system (12.10)) during the time evolution. If they are desired variables for 
output, one only needs to convert between p, e and z, T at the final output time. 

4.1. Asymptotic Property of the New Scheme. In this subsection we show that the new scheme, 
when applied to the quantum BGK equation, has the property (| 1 . 6[) . Consider the following time 
discretization: 

(46) r +1 - r l v Vxfn - n - km- n | \(m- +i 

At e e 

Some simple mathematical manipulation on (|4.6I) gives 
(4.7) 

1 + (X — 1)— At A— 

f n + l_ M n+l = y )e {f n_ M n ) _^^ v , VJ n + {M n_ M n + l ) + _^_ {M n + l_ M n^ 



Assume all the functions are smooth. When A > |, 



(4.8) \f n+1 -M n q +1 \ <a\f n -M%\ + 0(e + At), 

where < a — |1 + (A — 1) — 1/|1 + A— | < 1 uniformly in e and At. The 0(e) term comes from the 
second term of the right hand side of (|4.7I) . The 0{At) term is from the third and fourth terms. Then 

(4.9) \f H - M n q \ < a n \f° - M° q \ + 0(c + At) . 

Since At is taken bigger than e, this implies the property (11.61) . It is interesting to point out that / 
approaches M q , not M c , with (|4.6p . 

Remark 4.1. The first order (in-time) method can be extended to a second order by an Implicit- 
Explicit (IMEX) method (see also [B]^: 



(4.10) 



r-r „ fn _ gjCT) - KM n c - n \(M* C - n 

At/2 +V - V *t - e + 

f n+1 -r , v >v = Q q (n - hm* c - n , x(M n c - n + x(M n c +i - r +i ) 



At J e 2e 

This scheme can be shown to have the same property 11.6}) on the quantum BGK equation 
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5. Numerical Examples 



In this section, we present some numerical results of our new scheme (|4.5p (a second order finite 
volume method with slope limiters [8] is applied to the transport part) on the 1-D shock tube problem. 
The initial condition is 



The particles are again assumed to be the 2-D Maxwellian molecules and we adjust 9$ to get 
different initial data for both the Bose gas and the Fermi gas. 

In all the regimes, besides the directly computed macroscopic quantities, we will show the fugacity 
z and temperature T as well. They are computed as follows. First, (|2.10p (d v — 2) leads to 



We treat the left hand side of (|5.2[) as one function of z, and invert it by the secant method. Once z 
is obtained, T can be computed easily using for example the first equation of (|2.10[) . To evaluate the 
quantum function Q u (z), the expansion (|2 . 13[) is used for the Bose-Einstein function. The Fermi-Dirac 
function is computed by a direct numerical integration. The approach adopted here is taken from |15] 
(Chapter 6.10). 

When approximating the collision operator Q q , we always take M = 4, TV = 32 and L = 8, except 
L — 6 for the Bose gas in the quantum regime. 

5.1. Hydrodynamic Regime. We compare the results of our new scheme (I4.5[) with the kinetic 
scheme (KFVS scheme in [7]) for the quantum Euler equations (|2.9I) . The time step At is chosen by 
the CFL condition, independent of e. FigJT] shows the behaviors of a Bose gas when 9q = 0.01. Figj8] 
shows the behaviors of a Bose gas when 6q — 9. The solutions of a Fermi gas at #o = 0.01 are very 
similar to FigH so we omit them here. Fig|§] shows the behaviors of a Fermi gas when 8q = 9. All 
the results agree well in this regime, which exactly implies the scheme (|4.5p is asymptotic preserving 
(when the Knudscn number e goes to zero, the scheme becomes a fluid solver). 

5.2. Kinetic Regime. We compare the results of our new scheme (|4.5[) with the explicit forward 
Euler scheme. The time step At for the new scheme is still chosen by the CFL condition. When the 
Knudsen number e is not very small, lO^ 1 or 10~ 2 , the above At is also enough for the explicit scheme. 
FigfTUl shows the behaviors of a Bose gas when 9o = 0.01. FigfTTI shows the behaviors of a Bose gas 
when 6*o = 9. The solutions of a Fermi gas at 0$ = 0.01 are very similar to FigfTUl so we omit them 
here. FigfT2l shows the behaviors of a Fermi gas when 9q = 9. Again all the results agree well which 
means the scheme (|4.5[) is also reliable in the kinetic regime. To avoid the boundary effect, all the 
simulations in this subsection were carried out on a slightly larger spatial domain x £ [—0.25, 1.25]. 



A novel scheme was introduced for the quantum Boltzmann equation, starting from the scheme in 
The new idea here is to penalize the quantum collision operator by a 'classical' BGK operator so 
as to avoid the difficulty of inverting the nonlinear system p — p{z, T), e = e(z, T). The new scheme is 
uniformly stable in terms of the Knudsen number, and can capture the fluid (Euler) limit even if the 
small scale is not numerically resolved. We have also developed a spectral method for the quantum 
collision operator, following its classical counterpart |T2l 15]. 

So far we have not considered the quantum gas in the extreme case. For example, the Bose gas 
becomes degenerate when the fugacity z — 1. Many interesting phenomena happen in this regime. 
Our future work will focus on this aspect. 



(5.1) 




(5.2) 



Ql(z) = «oP 



6. Conclusion 
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Figure 7. Bose gas. e = le - 4, 6 = 0.01, z x = 0.0016, z r = 7.9546e - 04. Density 
p, velocity u, fugacity z and temperature T at t = 0.2. Af = 0.0013, Ax — 0.01. 
Solid line: KFVS scheme [7] for quantum Euler equations (|2.9p : o: New scheme (|4.5I) 
for QBE (TO) . 
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u, fugacity z and temperature T at t = 0.2. At = 0.0017, Ax — 0.01. Solid line: 
KFVS scheme [7] for quantum Euler equations (|2.9I) : o: New scheme (|4.5I) for QBE 
(ED). 
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Figure 9. Fermi gas. e = le - 4, 8 Q = 9, zi = 3.1887, z r = 1.0466. Density p, 
velocity u, fugacity z and temperature T at t = 0.2. At = 0.0013, Ax = 0.01. Solid 
line: KFVS scheme [7] for quantum Euler equations (12.91) ; o: New scheme (|4.5[) for 
QBE (fTTTT) . 
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Figure 10. Bose gas. e = le-2, 6 a = 0.01, zi = 0.0016, z r = 7.9546e- 04. Density 
p, velocity u, fugacity z and temperature T at t — 0.2. At — 0.0013, Ax = 0.01. 
Solid line: Forward Euler scheme for QBE (H~T|) ; o: New scheme (|4"3)l for QBE (Tl~T|) . 
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Figure 11. Bose gas. e = le - 1, 9 = 9, z t = 0.7613, z r = 0.5114. Density p, 
velocity u, fugacity z and temperature T at t = 0.2. At = 0.0017, Aa; = 0.01. Solid 
line: Forward Euler scheme for QBE (fTTj) ; o: New scheme (|4.5p for QBE (jl.ip . 
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Figure 12. Fermi gas. e = le - 2, O = 9, z\ = 3.1887, z r = 1.0466. Density p, 
velocity it, fugacity z and temperature T at t = 0.2. At = 0.0013, Aa; = 0.01. Solid 
line: Forward Euler scheme for QBE (jl.lj) : o: New scheme (|4.5J) for QBE (II. lj) . 
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